{r} [Zrodlo] danych] (http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx) {r} [Opis danych] (https://www.nature.com/articles/s42256-020-0180-7)
raport powinien zaczynać się od rozdziału podsumowującego całą analizę, streszczającego najważniejsze spostrzeżenia analityka
Kod wyliczający wykorzystane biblioteki. Kod zapewniający powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych. Kod pozwalający wczytać dane z pliku. Kod czyszczący dane (np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych). Sekcję podsumowującą rozmiar zbioru i podstawowe statystyki. Szczegółową analizę wartości atrybutów. Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji. Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie. Sekcję próbującą stworzyć klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji. Analizę ważności atrybutów najlepszego znalezionego modelu. Dodatkowo punktowane będzie wykonanie analizy typowej dla danych klinicznych, np. regresji logistycznej wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors) lub regresji Coxa (ang. Cox Proportional-Hazards Model).
Analityk nie musi, a nawet nie powinien, ograniczać się do powyższych punktów. Wszelkie dodatkowe techniki analizy danych, wizualizacje, spostrzeżenia będą pozytywnie wpływały na ocenę.
Ewentualne konkluzje, znalezione zależności warto potwierdzić dokonując sprawdzenia istniejących wyników badań w literaturze naukowej (np. na Google Scholar czy PubMed).
TBA
library(xlsx)
library(DT)
library(knitr)
library(dplyr)
library(tidyr)
library(janitor)
library(imputeTS)
library(lares)
library(plotly)
library(caret)
filename <- 'res/wuhan_blood_sample_data_Jan_Feb_2020.xlsx'
raw_data <- read.xlsx(filename, 1)
raw_data <- as_tibble(raw_data)
# OK
#how many NA in variables
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596
# var min i max
min(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] -1
max(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] 50000
glimpse(raw_data)
## Rows: 6,120
## Columns: 81
## $ PATIENT_ID <dbl> 1, NA...
## $ RE_DATE <dttm> 2020...
## $ age <dbl> 73, 7...
## $ gender <dbl> 1, 1,...
## $ Admission.time <dttm> 2020...
## $ Discharge.time <dttm> 2020...
## $ outcome <dbl> 0, 0,...
## $ Hypersensitive.cardiac.troponinI <dbl> NA, N...
## $ hemoglobin <dbl> NA, 1...
## $ Serum.chloride <dbl> NA, N...
## $ Prothrombin.time <dbl> NA, N...
## $ procalcitonin <dbl> NA, N...
## $ eosinophils... <dbl> NA, 0...
## $ Interleukin.2.receptor <dbl> NA, N...
## $ Alkaline.phosphatase <dbl> NA, N...
## $ albumin <dbl> NA, N...
## $ basophil... <dbl> NA, 0...
## $ Interleukin.10 <dbl> NA, N...
## $ Total.bilirubin <dbl> NA, N...
## $ Platelet.count <dbl> NA, 1...
## $ monocytes... <dbl> NA, 1...
## $ antithrombin <dbl> NA, N...
## $ Interleukin.8 <dbl> NA, N...
## $ indirect.bilirubin <dbl> NA, N...
## $ Red.blood.cell.distribution.width. <dbl> NA, 1...
## $ neutrophils... <dbl> NA, 6...
## $ total.protein <dbl> NA, N...
## $ Quantification.of.Treponema.pallidum.antibodies <dbl> NA, N...
## $ Prothrombin.activity <dbl> NA, N...
## $ HBsAg <dbl> NA, N...
## $ mean.corpuscular.volume <dbl> NA, 9...
## $ hematocrit <dbl> NA, 3...
## $ White.blood.cell.count <dbl> NA, 3...
## $ Tumor.necrosis.factorÎ. <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.concentration <dbl> NA, 3...
## $ fibrinogen <dbl> NA, N...
## $ Interleukin.1β <dbl> NA, N...
## $ Urea <dbl> NA, N...
## $ lymphocyte.count <dbl> NA, 0...
## $ PH.value <dbl> 7.415...
## $ Red.blood.cell.count <dbl> NA, 4...
## $ Eosinophil.count <dbl> NA, 0...
## $ Corrected.calcium <dbl> NA, N...
## $ Serum.potassium <dbl> NA, N...
## $ glucose <dbl> NA, N...
## $ neutrophils.count <dbl> NA, 2...
## $ Direct.bilirubin <dbl> NA, N...
## $ Mean.platelet.volume <dbl> NA, 1...
## $ ferritin <dbl> NA, N...
## $ RBC.distribution.width.SD <dbl> NA, 4...
## $ Thrombin.time <dbl> NA, N...
## $ X...lymphocyte <dbl> NA, 2...
## $ HCV.antibody.quantification <dbl> NA, N...
## $ D.D.dimer <dbl> NA, N...
## $ Total.cholesterol <dbl> NA, N...
## $ aspartate.aminotransferase <dbl> NA, N...
## $ Uric.acid <dbl> NA, N...
## $ HCO3. <dbl> NA, N...
## $ calcium <dbl> NA, N...
## $ Amino.terminal.brain.natriuretic.peptide.precursor.NT.proBNP. <dbl> NA, N...
## $ Lactate.dehydrogenase <dbl> NA, N...
## $ platelet.large.cell.ratio. <dbl> NA, 3...
## $ Interleukin.6 <dbl> NA, N...
## $ Fibrin.degradation.products <dbl> NA, N...
## $ monocytes.count <dbl> NA, 0...
## $ PLT.distribution.width <dbl> NA, 1...
## $ globulin <dbl> NA, N...
## $ γ.glutamyl.transpeptidase <dbl> NA, N...
## $ International.standard.ratio <dbl> NA, N...
## $ basophil.count... <dbl> NA, 0...
## $ X2019.nCoV.nucleic.acid.detection <dbl> NA, N...
## $ mean.corpuscular.hemoglobin. <dbl> NA, 3...
## $ Activation.of.partial.thromboplastin.time <dbl> NA, N...
## $ High.sensitivity.C.reactive.protein <dbl> NA, N...
## $ HIV.antibody.quantification <dbl> NA, N...
## $ serum.sodium <dbl> NA, N...
## $ thrombocytocrit <dbl> NA, 0...
## $ ESR <dbl> NA, N...
## $ glutamic.pyruvic.transaminase <dbl> NA, N...
## $ eGFR <dbl> NA, N...
## $ creatinine <dbl> NA, N...
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596
(np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych)
# NOT OK
#ogarnac te brakujace daty -- tyle ze ich nie ma?
# kolumny z samimy na
#replace -1 with NA
raw_data[raw_data==-1]<-NA
#PATIENT_ID
id_filled <- raw_data %>% fill(PATIENT_ID)
#OTHER NA VALUES
#is there any row that doesn't have NA value?
full_row_count <-dim(id_filled[complete.cases(id_filled),])[1]
print(full_row_count)
## [1] 0
# NA values in variables:
mean(is.na(id_filled %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770748
# there are two patients with negative days, as their only blood sample results arrived one day after their clinical outcome
#remove rows where all variables are empty
vars <- colnames(id_filled)[-(1:7)]
no_empty_rows<- id_filled[rowSums(is.na(id_filled[vars])) != length(vars), ]
no_empty_cols <- no_empty_rows[colSums(!is.na(no_empty_rows)) > 0]
# NA values in variables:
mean(is.na(no_empty_rows %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8658041
colnames_cleaned <- no_empty_cols %>% clean_names()
NA Values:
clean_NA<-function(column){
not_NA_count<-sum(!is.na(column))
if (not_NA_count>=2){
column <- na_interpolation(column, option = "linear")
column
}
else if (not_NA_count==1){
val <- first(na.omit(column))
column[is.na(column)] <- val
column
}
column
}
cleaned<- colnames_cleaned%>% group_by(patient_id) %>% mutate_each(funs(clean_NA))
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `mutate_each_()` is deprecated as of dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
var_df<-cleaned[-(1:7)]
“Sekcja podsumowująca rozmiar zbioru i podstawowe statystyki” Szczegółowa analiza wartości atrybutów.
DT::datatable(cleaned, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))
## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html
all_dim <-dim(cleaned)
tmp <- cleaned %>% select(patient_id, outcome, gender) %>% group_by(patient_id) %>% summarise(outcome_count = mean(outcome), gender_count = mean(gender))
## `summarise()` ungrouping output (override with `.groups` argument)
# ile pacjentow
patients_count <- length(distinct(cleaned, patient_id)$patient_id)
# ile pomiarow
measurements_count <- all_dim[1]
# ile srednio pomiarow na pacjenta
mean_measures <- measurements_count/patients_count
# ile kobiet/mezczyzn
genders <- tmp %>% count(gender_count) #Male 224, Female 151 # WYKRES
# ile umarlo/przezylo
outcomes <- tmp %>% count(outcome_count) #201 recovered, 174 died # WYKRES
# ile kolumn
columns_count <- all_dim[2]
# ile atrybutow
vars_count <- all_dim[2]-7
# ile zostalo wartosci NA
na_left <- round(100*mean(is.na(cleaned)), 0)
titles <- c("Liczba pacjentów", "Liczba pomiarów", "Średnia liczba pomiarów", "K", "M", "Smierc", "Zycie :(", "Liczba wierszy", "Liczba zmiennych", "Brakujace wartosci [%]")
values <- c(patients_count,measurements_count,mean_measures,1, 2, 0,1, columns_count,vars_count, na_left)
info_table <- data_frame(
titles,
values
)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
knitr::kable(info_table)
| titles | values |
|---|---|
| Liczba pacjentów | 361.00000 |
| Liczba pomiarów | 5606.00000 |
| Srednia liczba pomiarów | 15.52909 |
| K | 1.00000 |
| M | 2.00000 |
| Smierc | 0.00000 |
| Zycie :( | 1.00000 |
| Liczba wierszy | 80.00000 |
| Liczba zmiennych | 73.00000 |
| Brakujace wartosci [%] | 6.00000 |
cleaned %>% distr(outcome, gender)
#wykres death_bothgenders, alive_bothgenders
invisible(remove(tmp))
plot_data <- cleaned%>%select(patient_id, gender, admission_time, discharge_time, outcome) %>%group_by(patient_id) %>% summarise_all(funs(first))
admission_plot <- ggplot() +
coord_cartesian() +
scale_color_hue() +
facet_wrap(~outcome) +
layer(data=plot_data,
mapping=aes(
x=admission_time,
y=discharge_time, colour=factor(gender)
),
stat="identity",
geom="point",
position=
position_jitter()
)
ggplotly(admission_plot)
Szczegółową analizę wartości atrybutów.
#data.frame(unclass(summary(tmp)), check.names = FALSE, stringsAsFactors = FALSE)
res <- var_df%>% sapply(function(x){c(min=min(x,na.rm = TRUE), median=median(x,na.rm = TRUE), mean = mean(x,na.rm = TRUE),max=max(x,na.rm = TRUE), variance=var(x,na.rm = TRUE), standard_deviation=sd(x,na.rm = TRUE), na_count = sum(is.na(x)))}) %>% data.frame()
rd<- dim(res)
rows<-rd[1]
cols<-rd[2]
plot_data<-data.frame(c(1:cols),c(1:cols))
DT::datatable(res, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))
#plt<-ggplot(gather(var_df), aes(value)) +
# geom_histogram(bins = 10) +
# facet_wrap_paginate(~key,nrow = 1,
# ncol = 1,
# scales = "free")
#ggplotly(plt)
page_plot<-function(page_no, prow=3, pcol=3){
print(
ggplot(gather(var_df), na.action="na.omit", aes(value)) + ggtitle(" title")+
geom_histogram(bins = 10) +
facet_wrap_paginate
(~key, ncol = pcol, nrow = prow, scales = "free_x", page = page_no))
}
plot_row = 4
plot_col = 4
pages<- ceiling(cols/(plot_row*plot_col))
for (i in 1:pages){
page_plot(i, plot_row, plot_col)
}
“Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji.”
# OK
#TODO wiecej wynikow
plot_data <- corr_cross(var_df,
max_pvalue = 0.05,
top = 100,
rm.na=TRUE)
## Returning only the top 100. You may override with the 'top' argument
#plot(plot_data)
ggplotly(plot_data)
Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie.
timeline_plot <- ggplot() +
coord_cartesian() +
scale_color_hue() +
layer(data=cleaned,
mapping=aes(
x=re_date,
y=hypersensitive_cardiac_troponin_i, group=patient_id
),
stat="identity",
geom="point",
position=
position_jitter()
)
ggplotly(timeline_plot)
timeline_plot <- ggplot(cleaned, aes(x=re_date, y=serum_chloride, colour=factor(patient_id), group=patient_id)) + geom_line() + geom_point() + facet_wrap(~outcome)
ggplotly(timeline_plot)
#timeline_data <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose)
timeline_data <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose, outcome, gender) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose, outcome, gender) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose), outcome=first(outcome), gender=first(gender)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose, outcome, gender)
## `summarise()` regrouping output by 'id' (override with `.groups` argument)
timeline_plot <- ggplot(timeline_data,aes(day_no,avg_hemoglobin, shape=factor(gender)))+geom_point(color="blue")+geom_point(aes(day_no,avg_glucose, shape=factor(gender)),color="black") + facet_wrap(~outcome)
ggplotly(timeline_plot)
## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.
klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji.
ml_data<- cleaned%>%group_by(patient_id)%>%summarise_all(funs(last))
ml_data<-na_mean(ml_data) #TO ZROBIC OSOBNO DLA ZBIORU TESTUJACEGO I UCZACEGO!
#rozwazyc: usuniecie tych kolumn (wierszy), w których jest dużo wartosci NA (np. powyżej 40%?)
#ml_data$outcome=as.factor(ml_data$outcome)
ml_data$outcome=factor(ml_data$outcome,
labels = make.names(c("negative", "positive")))
set.seed(23)
inTraining <-
createDataPartition(
# atrybut do stratyfikacji
y = ml_data$outcome,
# procent w zbiorze uczącym
p = .75,
# chcemy indeksy a nie listę
list = FALSE)
training <- ml_data[ inTraining,]
testing <- ml_data[ -inTraining,]
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fit <- train(outcome ~ .,
data = training,
method = "rf",
trControl = ctrl,
ntree = 10)
rfClasses <- predict(fit, newdata = testing)
confusionMatrix(data = rfClasses, testing$outcome)
## Confusion Matrix and Statistics
##
## Reference
## Prediction negative positive
## negative 48 0
## positive 0 41
##
## Accuracy : 1
## 95% CI : (0.9594, 1)
## No Information Rate : 0.5393
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.5393
## Detection Rate : 0.5393
## Detection Prevalence : 0.5393
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : negative
##
ml_data$outcome=factor(ml_data$outcome,
labels = make.names(levels(ml_data$outcome)))
rfGrid <- expand.grid(mtry = 10:30)
gridCtrl <- trainControl(
method = "repeatedcv",
summaryFunction = twoClassSummary,
classProbs = TRUE,
number = 2,
repeats = 5)
set.seed(23)
fitTune <- train(outcome ~ .,
data = training,
method = "rf",
metric = "ROC",
preProc = c("center", "scale"),
trControl = gridCtrl,
tuneGrid = rfGrid,
ntree = 30)
rfTuneClasses <- predict(fitTune,
newdata = testing)
confusionMatrix(data = rfTuneClasses,
testing$outcome)
## Confusion Matrix and Statistics
##
## Reference
## Prediction negative positive
## negative 48 0
## positive 0 41
##
## Accuracy : 1
## 95% CI : (0.9594, 1)
## No Information Rate : 0.5393
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.5393
## Detection Rate : 0.5393
## Detection Prevalence : 0.5393
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : negative
##
ggplot(fitTune) + theme_bw()
### Analizę ważności atrybutów najlepszego znalezionego modelu
analiza typowa dla danych klinicznych, np.: